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ABSTRACT 

In  this  report  we  review  and  present  a  performance  comparison  of  a  number  of  recenfly 
proposed  range-based  source  localization  algorithms  that  approximate  the  maximum 
likelihood  estimator.  Although  most  of  these  algorithms  are  available  in  the  literature  we  give 
an  elaborated  presentation  of  fhem  and  provide  some  comments  on  their  performances.  The 
optimization  fechniques  used  in  fhese  localizafion  mefhods  range  from  convex  relaxation, 
global  continuation  to  the  generalized  trust  region  subproblem  method.  Among  the 
localization  algorithms  considered,  the  range  weighted  squared  range  least  squares  estimator 
(RW-SRLS)  seems  to  be  a  particularly  attractive  algorithm  as  it  is  fast  and  yet  only  marginally 
suboptimal  in  a  maximum  likelihood  sense,  especially  for  small  range  measuremenf  errors. 
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Recent  Advances  in  Source  Localization  Using  Range 

Measurements 


Executive  Summary 


In  this  report  we  review  and  present  a  performance  comparison  of  a  number  of 
recenfly  proposed  range  or  time  of  arrival  (TOA)-based  source  localization  algorifhms 
fhaf  approximafe  fhe  maximum  likelihood  estimator.  Some  of  fhese  algorifhms  are 
available  in  fhe  literature  and  others  are  contributed  for  fhe  firsf  time  in  this  report. 
This  comparison  study  is  important  because  TOA-based  geolocation  is  gaining 
popularity,  mainly  because  of  the  steady  emergence  of  more  stable  and  accurate  time 
clocks,  which  enable  dispersed  platforms/ vehicles/  soldiers  to  insfanfly  deduce 
separation  disfances  from  signal  flighf  times.  A  unif  is  fhen  able  to  rapidly  seK-localize 
in  difficult  environments  using  as  little  as  two  or  three  time/ position  stamped  signals 
(i.e.  transmit  time  and  position  information  is  included  in  the  signal's  preamble). 
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1.  Introduction 

The  passive  localization  of  emitters  is  of  considerable  interest  in  many  applications 
including  defence,  wireless  communications  and  navigation.  A  number  of 
localization  techniques  were  proposed  in  the  past  depending  on  the  type  of 
measurements  and  position  estimation  approach.  Most  popular  among  these  are 
the  time  of  arrival  (TOA),  time  difference  of  arrival  (TDOA)  and  angle  of  arrival 
(AO A)  localization  algorithms  (see,  e.g.,  [1,2,3,4,5,6]). 

In  this  paper  we  focus  on  range-based  localization  techniques  in  the  plane  which 
underpin  TOA  localization  and  are  extensively  utilised  in  sensor  network 
localization.  In  particular  we  compare  a  number  of  advanced  source  location 
estimators  and  examine  how  closely  they  approximate  the  least  squares  estimator. 
The  difference  between  these  estimators  is  either  the  objective  function  or  the 
computation  method  used  to  minimise  it.  The  two  most  popular  cost  functions 
considered  in  the  literature  are  the  range  least  squares  (R-LS)  and  the  squared 
range  least  squares  (SR-LS)  [3].  The  R-LS-based  formulation  is  of  great  interest  and 
has  been  known  for  its  optimal  performance  [3]  [7].  If  the  measurement  errors  are 
zero-mean  Gaussian,  then  it  coincides  with  the  maximum  likelihood  estimator 
(MLE).  The  main  challenge  is  how  to  efficiently  compute  an  R-LS  position  estimate. 

A  number  of  optimization  tools  may  be  applied  to  globally  solve  the  R-LS  problem 
and  are  usually  computationally  intensive.  Polynomial  continuation  methods  [8]  or 
computer  algebra  techniques  [9]  [10]  can  determine  the  global  minimiser  after 
transforming  the  R-LS  problem  into  a  system  of  multivariate  polynomial 
equations.  Other  solution  methods  are  based  on  semidefinite  relaxation  (SDR) 
[7]  [11]  [12]  [13]  or  polynomial  optimization  (POP)  [14].  They  can  approximate  the  R- 
LS  estimator  to  an  arbitrary  degree  but  at  the  expense  of  a  significant  increase  in 
computer  processing  and  memory  storage.  Another  approach  to  address  the  R-LS 
problem  is  to  apply  a  smoothing  kernel  to  the  LS  cost  function  and  progressively 
minimise  it  by  continuation  [15]  [16]. 

In  the  first  part  of  the  paper  we  present  source  localization  approaches  based  on 
squared  range  measurements.  The  SR-LS  source  position  estimator  analysed  in  [3], 
is  one  classic  estimator,  which  exhibits  several  nice  mathematical  properties.  The 
position  estimate  can  be  computed  exactly  and  efficiently  using  generalized  trust 
region  subproblem  (GTRS)  optimization  techniques  [17]  [18]  [19].  However,  as  we 
shall  see  later  in  the  paper,  the  estimator  is  not  optimal  in  the  maximum  likelihood 
sense  [3].  We  propose  a  variant  of  this  estimator  which  uses  range-based  corrective 
weights  to  reduce  the  position  estimation  bias  error.  Most  importantly  the 
computational  efficiency  of  this  variant  estimator,  which  we  call  range  weighted 
SR-LS  (RW-SRLS),  is  the  same  as  that  of  the  SR-LS  [3].  Computer  simulations  of  the 
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RW-SRLS  estimator  show  definite  advantages  in  terms  of  position  accuracy  and 
low  computational  effort. 

In  the  second  part  of  this  paper  we  briefly  present  and  compare  the  performance  of 
some  of  the  algorithms  known  to  solve  or  best  approximate  the  R-LS  estimator. 
These  include  two  contributions  based  on  computer  algebra  and  convex  relaxation. 
The  earliest  research  effort  to  mathematically  approximate  the  R-LS  solution  was 
made  by  Cheung  et  al.  [7].  Their  proposed  approximation  method  is  a 
straightforward  convex  relaxation  of  the  R-LS  problem.  Further  work  in  [3] 
provided  a  rigorous  analysis  of  their  algorithm.  As  demonstrated  in  [3]  its 
performance  is  often  degraded  because  the  relaxation  solution  can  be  significantly 
different  from  that  of  the  original  R-LS  problem.  An  alternative  relaxation 
formulation,  called  CPnBSP-SDR,  is  proposed  in  this  paper,  which  shows 
improved  performance  accuracy.  This  work  is  related  to  that  of  [13]  and  both  are 
computationally  more  efficient  than  the  SLCP  positioning  method  [12],  which  is 
restricted  to  two-dimensional  geolocation  problems.  We  next  investigate  the  kernel 
continuation  method  of  [16],  which  offers  some  performance  compromise  by  using 
a  smoothing  parameter  to  control  the  shape  of  the  cost  function  and  progressively 
steer  the  solution  towards  the  global  minimum.  However,  it  remains  unclear  how 
to  update  the  smoothing  parameter  and  a  global  convergence  certificate  is  not 
provided.  Furthermore  the  continuation  derivations  in  [16]  are  only  applicable  to 
source  localization  problems  in  the  plane.  Finally  we  provide  an  algorithm  which 
solves  the  R-LS  problem  exactly  subject  only  to  numerical  round-off  errors.  This 
algorithm  is  based  on  the  work  [10]  and  transforms  the  problem  of  solving  a 
multivariate  polynomial  system  into  the  null  space  analysis  of  a  large  sparse 
matrix  whose  nonzero  entries  consist  of  shifted  polynomial  coefficients.  This 
method,  however,  is  suitable  for  smaller  localization  problems  in  the  plane  as  it 
becomes  computationally  prohibitive  when  the  number  of  measurements  increases 
beyond  three. 

This  report  comprises  a  number  of  sections.  Section  2  introduces  the  R-LS  problem 
and  provides  some  basic  results  related  to  the  CRLB  and  the  regional  bound  for  the 
R-LS  solution.  The  rest  of  the  report  is  mainly  divided  into  two  parts  where  the 
first  covers  two  variants  of  the  SR-LS  approach  (Section  3).  The  second  part 
(Sections  4  and  5)  covers  R-LS  relaxation  and  continuation  based  source 
localization  approaches.  Section  6  presents  Reid  and  Zhi's  null  space  method  to 
solve  the  R-LS  problem  exactly.  Section  7  presents  several  important  results 
highlighting  performance  comparisons  among  various  source  position  estimators. 
Finally  concluding  remarks  are  made  in  Section  8. 
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2.  Problem  Formulation  and  Basics 

2.1  Problem  Formixlarion 


In  tills  subsection  we  introduce  tlie  source  localization  problem  to  be  solved. 
Unless  specifically  mentioned  we  express  tlie  somce  localization  problem  in  91", 
where  » =  2  and  ??  =  3  refer  fo  source  geolocafion  in  die  plane  and  in  space, 
respectively.  As  will  become  cleai'  later  on,  some  localization  algoiitlmis  are 
applicable  to  both  91“  and  91  \  wliile  odiers  are  only  apphcable  in  9l’.  However, 
all  computer  simulation  results  are  restiicted  to  91^.  If  we  assume  diaf  the  range 
measurement  eiTors  are  independent  zero-mean  Gaussian,  then  the  maximrmi 
likelihood  (ML)  somce  position  estimator  reduces  to  die  noidineai’  LS  problem, 
wliich  is  denoted  bv  R-LS  as  given  in  [3] 


P 


Figure  1:  Emittei-receivei  geoinetnj.  The  eniittei  is  at  an  unknown  position,  p,  to  be  estimated 
given  measurements  of  the  distances  separating  it  from  the  sensors. 


uiin  /(P)  =  X’''.(||P-P,||-^,)'  (1) 

As  dlustrated  in  Figure  1,  the  points,  p, ,  are  die  known  positions  of  N  sensors  in 
91"  and  d,  their  associated  measmed  ranges  to  the  somce  location.  Tlie  non- 

N 

negative  numbers,  w. ,  aie  positive  normalized  weights  (i.e.  =  1 ).  Tlie  weights 

1=1 

are  usually  set  proportional  to  the  reciprocal  of  respective  range  measurement 
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error  variances.  For  a  localization  problem  in  the  plane  (i.e.  n  =  2),  the  coordinates 
of  the  source  position,  p,  are  denoted  as  (x,  They  become  {x,y,zy  for  a 
localization  problem  in  space  (i.e.  n  =  3). 

Before  reviewing  a  number  of  range-based  source  localization  works,  we  first 
establish  a  bound  for  the  solution  of  (1). 

2.2  Basic  Results 


The  gradient  at  a  minimum  of  (1)  vanishes  leading  to 

Vp/(p)  =  X’^/t-P/ (2) 

/=1 

where  u,  =  jp — ^  is  a  unit  direction  vector  pointing  from  the  i*  sensor  position  to 

Ip-p-II 

p .  Expanding  (2)  leads  to 

p  =  Z^-(p^+^“J  (3) 

/=1 

N 

If  we  let  Pc  =  X  ^'Pi '  sensor  position  centroid,  then  (3)  reduces  to 

/=1 

N  ^ 

p  =  p^  +  ^iv.d,u,-  (4) 

z=l 

The  next  theorem  shows  that  p  is  confined  within  the  ball, 

={pe9t",||p-Pc||<i?}  (5) 

N  ^ 

where  R  -  . 

/=! 

Theorem  1:  An  optimal  solution,  p,  of  (1)  is  bounded  within  a  ball  centred  at 

N  ^  ^ 

and  of  radius  R-^  w^.  . 

i=l  i=l 


Proof: 


Let  Pc  =  Z^;P;  the  sensor  centroid  position.  From  (4)  an  optimal  solution  to  (1) 

!=1 

N  ^ 

must  satisfy  p  -Pc  =  ^PPty  Euclidean  norm  to  both  sides  of  this 


equation,  we  obtain  p  -  Pc  = 


Zw.d,u. 


Z=1 


resulting  in  ||p  -  Pc  ||  <  £  IK'  II  =  Z  • 


This  means  that  p  belongs  to  the  ball  ^  =  |p  e  W,  p  -Pc  ^  r}- 


i=l 

□ 
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Finally  the  error  performance  of  all  geolocation  estimators  needs  to  be  assessed 
against  the  Cramer  Rao  Lower  Bound  (CRLB).  Under  the  Gaussian  range 
measurement  error  assumption,  R-LS  is  an  asymptotically  efficient  and  unbiased 
estimator  as  N  ^  co,  achieving  an  estimation  covariance  identical  to  the  Cramer- 
Rao  lower  bound  (CRLB): 


CRLB=.|  |;^(p-p,)(p-p,)"  (6) 

\  1=1  cr,-  dj 

where  d-  -  lip  -p,  ||  and  a.  is  the  i*  range  noise  standard  deviation.  Even  though  R- 


LS  is  considered  to  be  a  benchmark  because  of  its  asymptotic  efficiency  [4],  it  is 
neither  efficient  nor  unbiased  for  finite  N. 


3.  Squared  Range  Based  Estimation  Approaches 

In  this  section  we  focus  on  source  localization  approaches  based  on  squared 
ranges.  The  main  difficulty  of  finding  the  LS  global  minimiser  is  the  presence  of  the 
square  root  operator  in  the  objective  function  (1).  If  we  square  the  norm  in  each 
term  and  compare  it  to  the  square  of  its  associated  measured  range,  then  a 
polynomial  objective  function  of  degree  four  emerges.  A  rich  repertoire  of 
polynomial  and  quadratic  optimization  techniques  may  then  be  applied  to  globally 
solve  the  problem.  However,  the  source  position  estimate  is  only  suboptimal  in  the 
maximum  likelihood  sense.  In  the  following  we  present  and  compare  two  squared 
range  based  least  squares  (SR-LS)  estimation  techniques. 

The  basic  SR-LS  source  position  estimate  is  obtained  by  solving 

ming(p)  =  ^w,(||p-p,f -JM  (7) 

where  the  parameters,  w.,  are  positive  weights.  In  what  follows  we  will  derive  a 
variation  of  this  estimator,  which  proves  to  be  more  optimal  in  a  maximum 
likelihood  sense 

3.1  Range  Weighted  SR-LS 


Observe  that  g(p)  in  (7)  can  also  be  expressed  as 

N  __  2  ^2 

g’(P)  -  (l|P“Pj“^i)  (l|P“Pi|l'''^>)  and  note  that  when  p  is  in  the  vicinity  of  the 

i=i 

estimated  source  position,  we  have  ||p-pjsid.  and  therefore 
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^  ^  2  /  -  \2 

g^(P)~2^^^/(||P“P,||“^,)  (^^i)  •  Clearly  fomiulation  (7)  weiglis  longer 
/=1 

measmement  ranges  significantly  more  tlian  shorter  ones.  Cancelling  tliis 
imdesirable  effect  leads  to  the  balanced  fomuilation 


wliich  is  more  optimal  in  a  maxiuuim  likehliood  sense.  We  denote  position 
estin\ator  (8)  as  RW-SRLS,  wliich  stands  for  Range  Weighted  SR-LS.  As  we  sliall 
see  later  in  the  computer  simulations  section  (Section  7),  the  RW-SRLS  estimator 
results  in  significant  improvement  in  localization  accmacy. 


Figure  2:  Contour  ylot  of  the  weighted  SR-LS  objective  function.  The  global  minimum  is  shown  ns  a 
red  diamond.  The  blue  t  inngle  coiresponds  to  a  local  minimum.  The  black  square  is  the 
tine  source  position. 
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Figure  3:  Contour  plot  of  tlie  R-LS  objective  function.  The  global  niininium  is  slmvn  as  a  red 
diamond.  The  blue  tiiangle  conesponds  to  a  local  miniinnm.  The  black  square  is  the  true 
source  position. 

Both  tlie  SR-LS  and  RW-SRLS  source  position  estimation  problems  are  instances  of 
the  generalized  trust  region  subproblems  (GTRS)  [17][18][19].  Tlrey  are 
tiansfoniied  to  quadratic  optimization  problems  sirbject  to  a  single  quadratic 
coirstrarnt  before  bemg  solved.  A  rigorous  application  of  tliis  approach  is  given  in 
[3],  wlrich  includes  formulation  (7)  but  not  (8).  For  more  details  orr  how  to  solve  (7) 
or  (8),  refer  to  [3]  [19]. 

Althouglr  (8)  is  a  closer  approximation  to  (1),  the  R-LS  arrd  RW-SRLS  global 
murimisers  may  occasionally  differ  considerably.  Tire  followirrg  nrmrerical  example 
iUrrstrates  tlris  point. 

Numerical  Example:  Consider  a  soruce  located  irr  the  plarre  at  positiorr  (11.24,  - 
5.32)  and  asstmre  the  availability  of  5  sensors  whose  locatiorr  coordinates  are  (7.37, 
4.36),  (11.00,  4.22),  (-8.89,  -6.49),  (10.20,  7.39)  and  (-9.84,  -6.87).  Let  tire  rroisy 
measmerrrerrts  be  r^espectively,  8.45,  11.67,  17.10,  13.32  arrd  19.81.  Hrese  are 
obtaiired  by  addiirg  zero  mearr  Gaussian  rroise  of  starrdar  d  deviatiorr  (cr  =  1.5)  to 
dre  true  rarrge  measuremerrts.  Figmes  2  arrd  3  show  the  corrtom’  plots  associated 
with  the  weighted  SR-LS  arrd  R-LS  objective  fmrctiorrs,  respectively.  Tire  weiglrted 
SR-LS  global  mirrinrrmr  iir  Figure  2,  showtr  as  a  red  diarrrorrd,  occms  at  (-2.97, 
10.19)  arrd  Iras  a  value  of  10.08.  Tire  figrrre  also  depicts  a  local  nrirrinrum  (showrr  as 
a  blue  triarrgle)  at  (11.52,  -5.08)  with  arr  objective  value  of  11.20.  Figme  3,  orr  the 
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other  hand,  displays  an  LS  global  minimum  of  10.95  at  (11.67,  -5.53).  The  local 
minimum  occurs  at  (-3.29,  10.35)  and  has  an  objective  value  of  11.21.  The  two 
figures  demonstrate  that  the  R-LS  global  minimiser  may  be  very  distant  from  the 
WR-SRLS  global  minimiser.  In  this  example  the  R-LS  global  minimum  occurs  in 
the  vicinity  of  the  true  source  position  but  this  is  not  always  true  in  general. 


4.  Geolocation  Using  Semidefinite  Relaxation  (SDR) 


In  the  following  we  present  three  R-LS  approximation  algorithms,  which  are  based 
on  the  semidefinite  relaxation  of  (1)  using  different  formulation  techniques. 
Solving  the  relaxed  problem  is  straightforward  but  the  obtained  solution  may  or 
may  not  correspond  to  that  of  (1).  Fortunately  a  convergence  certificate  to  the 
global  solution  of  (1)  is  provided,  in  the  form  of  a  matrix  rank  test. 

4.1  Cheung's  Method  [7] 


This  is  perhaps  the  first  work  which  applied  convex  relaxation  techniques  to  (1). 
The  algorithm  details  are  given  in  [7]  for  localization  problems  in  9?^ .  However, 
for  more  generality  and  completeness  we  briefly  reproduce  this  algorithm  using 
notations  applicable  to  both  9?^  and  9?^ . 


Problem  (1)  can  be  equivalently  re-expressed  as 


~ 


min  /(p)  =  £^-{5“<) 

P,-;  ^  V  / 


ri  =  P-P, 


l<i<N 


(9) 


If  we  define  G and  P  =  [p^,i]^[p^,i]g  then 

(9)  may  also  be  formulated  equivalently  as 


min  Ew,(G,-2G,,,d,+d,^) 

;=i 

G,=ir{C,P},  l</<iV 
[G  0l 

>0 

0  P 

p  =  p  =1 

^N+\,N+\  ^  n+\,n+\  ^ 

rank(G)  =  rank(P)  =  1 
where 


(10) 
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By  dropping  the  rank-1  constraint,  problem  (10)  is  relaxed  to 

min  ^  w.  (g,,  -  2 G^^i,  d,  +  d," ) 

;=1 

G,  =ir{C,P},  l</<iV 


[0  PJ 

^  jV+l,N+l  ^tt+1,/2+1  ^ 


(11) 


which  can  be  solved  efficiently  using  interior  point  methods.  It  should  be  pointed 
out  that  the  global  solution  of  (11)  coincides  with  that  of  (10)  only  if  the  matrices  G 
and  P  are  rank-1  (or  nearly  rank-1  in  practice).  The  solution  is  extracted  from  the 
last  column  (or  row)  of  P.  A  recent  result  due  to  [3]  demonstrated  that  the 
solution  G  of  the  relaxed  problem  is  guaranteed  to  be  rank-1.  This  property, 
however,  does  not  always  carry  over  to  P . 

When  the  solution  P  is  of  rank  larger  than  1,  which  is  not  rare  in  practice,  the 
global  solution  to  (10)  may  not  be  known,  but  the  computed  global  minimum 
remains  a  lower  bound  to  the  objective  function  of  (1),  (9)  or  (10).  In  practice  if  P  is 
not  of  rank  1,  a  coarse  approximate  solution  to  (1)  may  be  obtained  by  applying  the 
singular  value  decomposition  to  P  to  find  a  reduced  rank  approximation.  With 
more  effort  one  may  obtain  'an  improved'  solution  through  a  randomization 
process  [20].  This  is  briefly  implemented  as  follows.  We  first  exclude  the  last 
column  and  row  of  P  and  call  the  resulting  submatrix,  P^.  We  also  let  t  be  the  last 

column  of  P ,  without  the  last  entry,  '1'  .  The  matrix,  E  =  P^  -tt^,  is  an  nxn  error 
covariance  matrix.  One  may  randomly  pick  v  as  a  Gaussian  random  vector  from 
the  distribution  v  ~  A^(t,P^  -tt^),  because  v  solves  the  optimization  problem  (9)  on 
average.  If  enough  samples  are  taken  using  this  distribution,  then  one  may  obtain 
'an  improved'  solution  for  the  problem  by  selecting  the  sample  resulting  in  the 
lowest  objective  value  of  (1).  Alternatively  one  may  be  interested  in  finding  a 
region  where  the  source  lies  with  a  given  probability,  p.  The  covariance  matrix  E 
defines  an  elliptic  region  given  by  (v-t)^E“'(v-t)  <  d  ,  where  d  controls  the  size 
of  the  elliptic  region  as  a  function  of  p  .  For  a  localization  problem  in  the  plane,  the 
size  parameter,  d ,  is  related  to  the  probability,  p ,  through  S  --2  log(l  - p)  [21]. 

4,2  The  SLCP  algorithm  [12] 

The  source  localization  in  the  complex  plane  (SLCP)  is  another  semidefinite 
relaxation  program  applied  to  the  LS  cost  function  of  (1).  This  program  is  only 
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applicable  to  localization  problems  in  the  plane  (i.e.,  n  =  2).  Its  main  approach  is  to 
transform  the  localization  problem  from  9?^  to  the  complex  plane.  Any  point  on  a 
circle  can  be  represented  in  the  complex  domain  as  s,  =p,  ,  where  p,  is  the 

circle's  centre  and  its  radius.  The  i*  point  position  on  the  circle  is  indicated  by 
the  phase  angle,  . 


A  ^2 

Minimizing  the  cost  function  /(p)  =  ^  tv,  (||p  - p,  ||  -  d^  ^  is  equivalent  to  finding  the 


source  position  closest  to  the  circles  in  a  square  sense. 


mm 

P.S, 


in  Jw.||p-s,. 


(12) 


subject  to  ||p;-S;||  =  d; 

N 

It  can  be  shown  [12]  that  the  optimal  solution  to  (12)  satisfies  p*  =  ¥  =S 

z=l 

where  ¥  is  the  weighted  average  of  the  optimal  circle  points,  s* . 


In  the  complex  domain  problem  (12)  can  be  re-expressed  as 


N 

min  V  w,  \\s-s. 
^  dl 


subject  to  s.  =  p,  +  d^e^*' 

If  we  define  0  =  f  e  ,  2i  =  {p^,p2,.. 

D  =  diag{[d^,d2,--- ,df^f) ,  w  =  [^1,^2,...,  and  H  -  diag(vf)  -  vfvf^ , 

becomes 


(13) 


,Pn 


e  C 


then  (13) 


min  (a  +  De)"i:(a  +  De) 
subject  to  IbI  =  1 


(14) 


Expanding  (14)  and  eliminating  constants  lead  to 

min  2Re{c"0}-0"M0 
0 

subject  to  |0|  =  1 

where  =  a^ID  and  M  =  Dww^D 


(15) 


Examining  the  objective  function  of  (15)  we  observe  that  if  we  substitute  0  by  0e^“ , 
only  the  linear  term,  2Re{c^0},  changes.  This  term  reaches  its  minimum  of  -2|c^0| 

when  c^0  becomes  real  and  negative.  Hence,  problem  (15)  reduces  to 
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max  2  c  0  +0"M0 

6  I  I 

subject  to  |0|  =  1 


or 


max  2^tr{cc"QQ^}  +  tr{M00"} 


(16) 


(17) 


subject  to  |0|  =  1 

where  the  phase  vector,  0 ,  is  related  to  that  of  (16)  by  a  common  phase  shift. 


By  putting  O  =  00^  and  following  standard  semidefinite  programming  relaxation 
techniques,  (17)  becomes 


max  i  +  ir{M<D} 
subject  to  O  >  0 


(18) 


4ir{cc"<D}  > 

which  can  be  solved  efficiently  using  interior  point  methods.  In  real  form  (18) 
becomes 


max  i  +  irlMX} 

x,\,t 


subject  to 


L-Y  XJ 

X..  =1 
11 

4ir{Re{cc"}X-Im{cc"}Y}  > 


(19) 


where  X  and  Y  (skew-symmetric)  correspond  to  the  real  and  imaginary  parts  of 
O.  The  solution  may  be  obtained  using  singular  value  decomposition  of  O  (see 
[12]  for  more  details).  It  should  coincide  with  that  of  (1)  if  the  rank  of  O  is  1. 

4.3  Closest  Point  to  N  Ball  Surfaces  Problem 

In  this  section  we  present  an  improved  convex  relaxation  formulation  of  the  R-LS 
problem  (1),  which  is  applicable  to  both  9?^  and  9?^ .  A  similar  work  has  also 
appeared  in  [13]  but  we  will  present  our  variant  convex  relaxation  of  the  problem 
and  indicate  the  differences  later  in  this  section. 


The  gradient  at  a  minimum  of  (1)  is  given  by  (2),  where  u  ^  is  a  unit 

Ip-P/I 

direction  vector  pointing  from  the  i*  sensor  position  to  the  source. 
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Figure  4:  Geometiic  interpretation  oftlie  first  order  optimality  condition  in  the  plane. 


If  we  substitute  s,  =  p,  +  .  iu  (2),  then  from  Figme  4  it  is  clear  that  each  residual 

eiTor,  p-s,,  is  radial  to  the  circle/ sphere  given  by  ||s, -p,||" -d/  =  0.  Hris  is  iirdeed 


f 


\ 


true  because  p  -  s,  =  p  -  p,  -  <y,u,  = 


1-^ - !- 

V  iP-P.iiy 


(p-p,).  Furdiermore  the  i'^'  temi  in 


the  objective  frmction  can  be  re-expressed  as  a  frmction  of  the  residual  eiior  as 


In  Hght  of  tlris  we  may  write 


P-P, 


=  l|P-s.- 


N 

min  2^  m' 

p,s, 

subject  to  ||s,-  -  p,-||'  = 


HiP-S/ 


(20) 


Problem  (20)  may  be  looked  at  as  die  problem  of  finding  a  point,  p  in  91",  closest 
in  a  weighted  LS  sense  to  N  ball  smfaces,  dB^_  ,  1  <  /  <  iV .  Hence  it  is  convenient 

to  refer  to  problem  (20)  as  the  closest  point  to  iVball  smiaces  problem  (CF^'BSP). 
Writing  dowrr  the  first  optimality  condition  of  (20)  widr  respect  to  p  results  in 
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N 

which  shows  that  the  closest  point  to  A^ball  surfaces  is  a  weighted 

k=\ 

centroid  point  of  A^ball  surface  points  satisfying 

N  N 

minyw,  -s, 

/=1  k=\ 

subject  to  s.  ^  p .  +  d-u.,  (21) 

||u,f  =1,  l</<iV 

Clearly  both  problems  (20)  and  (21)  are  equivalent  and  the  only  difference  is  that  p 
is  made  explicit  in  (20).  A  similar  expression  to  (21)  appeared  in  [13].  The  unknown 

N  N  ^ 

variables  in  (21)  are  the  unit  vectors  u,.  Expanding  Xw,  results  in 

1=1  k=\ 

N  N  N 

Xw,.sfs; S  =  [si,s2,---,s^f,  w  =  [wi,W2,---,W/yf  and 

1=1  f=l  j=l 

L  =  ww^  -  J/ag(w)  then  this  latter  expression  can  be  written  compactly  as 
-ir|s^Ls|  where  tr{.}  refers  to  the  trace  operator.  If  we  further  let  a  =  [pl,p2,■■■,p^^f , 
D  =  diag([dj,d2,---,d^]  )  and  U  =  [uj,Uj,---,u^J  then  from  the  definition  of 
s .  =  p .  +  J,u .  it  follows  that  S  =  a  +  DU  and  (21)  becomes 

max  tr|  (a  +  DU)"'i:(a  +  DU)} 

(22) 

subject  to  |u|  =  1 

where  |U|  =  1  signifies  that  each  row  of  U  has  a  norm  of  1.  Note  that  although 

problem  (22)  appears  similar  to  (14),  formulation  (22)  is  not  restricted  to 
localization  problems  in  the  plane. 

Putting  c  =  DEa  and  M  =  Dww^D  and  expanding  (22)  leads  to  the  optimization 
problem 

max  tr\ c^U  +  U^c  +  U^Mul 

(23) 

subject  to  |u|  =  1 

Although  problem  (23)  is  difficult  to  solve  because  it  is  non-convex,  it  can  be 
shown  that  c^U  is  symmetric,  positive  semidefinite  at  an  optimal  solution.  We 
next  consider  relaxing  formulation  (23)  to  be  able  to  solve  it  or  find  an  approximate 
solution  to  it.  Let  T  be  the  augmented  n  -column  matrix  formed  by  appending  the 
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n X  n  identity  matrix,  ,  to  the  bottom  of  U  (i.e.  T  =  [U;I„ ] ).  If  we  foim  T  =  TT^, 

^  0 
I-. 

of  ^  is  7? .  Tlie  objective  fimction  of  (23)  may  be  expressed  in  tenns  of  ^  as 


then  T  can  be  par  titioned  as  ^  = 


where  =  UU^  and  0  =  U .  Tire  rank 


where 


M 


0. 


Usiirg  these  notations,  the  convex  relaxation  of  (23)  reads  as 


max 

^>0,  =1,  1<7<A^,  (24) 

V  =  I 

^  N+l:N+n.N+l:N+n 

where  the  non-convex  constraint,  lauk(^)  =  n  ,  Iras  beeir  dr  opped. 

Problem  (24)  carr  be  solved  efficierrtly  usirrg  irrterior  poiirt  methods.  Srnrilar  to  (23), 
it  carr  be  showrr  tlrat  c^0  is  positive  senridefurite  at  arr  optirrral  sohrtiorr.  Let  the 
solutiorr  of  (24)  be  'P* .  If  tire  rarrk  of  is  approxmrately  n  therr  0  is  rrumerically 
arr  approximate  solutiorr  of  (23).  Set  the  firral  sohrtiorr  to  0*  =  0 .  If  the  rarrk  is 
larger  tharr  n ,  therr  arr  approxinrate  somce  solutiorr  is  obtairred  as  follows 

1.  Compute  the  sirrgular  value  decompositiorr  (SVD)  of  as  L  A  .  Tire 
sirrgrrlar  values  iir  A  are  arrarrged  m  a  decreasirrg  order. 

2.  Extract  the  diagorral  77x77  submatrix.  A,  corrsistirrg  of  the  n  largest  sirrgular 
values  of  T*  (i.e.  the  n  largest  diagorral  values  of  A).  Sinrilariy  let  L  be  the 
77  cohmurs  of  L  associated  with  the  largest  n  smgular’  values,  after 
renrovirrg  the  bottom  n  rows.  Tire  size  of  L  is  x  77 .  Firrally  let  R  be  the  n 
colunurs  of  R  associated  with  the  largest  11  smgrrlar’  values  arrd  orrly 
keeprrrg  the  bottom  n  rows.  Tire  size  of  R  is  77  x  77 . 

3.  Compute  L  A  R^ ,  rrornralize  tire  rows  to  1  arrd  store  the  outcome  irr  0 . 

4.  Compute  the  SVD  of  A  =  c^0  as  K  and  retmir  0*  =  0QK^ . 

Tire  eshnrated  somce  locatiorr  follows  as 


p*=(a  +  D0*)^w  (25) 

Note  tlrat  irr  [13]  the  airthors  derived  a  slightly  differerrt  SDR  problem  to  (24).  Tire 
nrairr  differerrce  is  that  om'  relaxatiorr  formrrlatiorr  (24)  makes  explicit  use  of  0  arrd 
^  whereas  dreu  formulatiorr  (see  (16)  irr  [13]),  nrakes  use  of  ^  arrd  Z  where  Z  is 

related  to  VaA^  .  Tire  inrplicatiorr  of  tlris  is  drat  everr  if  the  rarrk  of  ^  is  n ,  the 
courprrted  U  (by  applyirrg  SVD  orr  ^> )  is  usually  rrot  the  problem  solutiorr  arrd  orre 
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has  to  adjust  U  tluough  application  of  an  optunal  orthogonal  traiisfomration  by 
executing  (12)  of  [13].  In  oxu'  work  tliis  step  is  needed  only  if  the  rank  of  is 
larger  than  n  and  is  caiiied  out  using  polar’  decomposition,  wliich  is  implemented 
using  SVD. 


5.  Geolocation  Using  the  Continuation  Method  [16] 


Tire  material  presented  in  tlris  section  is  mostly  based  orr  [16].  For  computation 
pruposes  we  deviate  horn  [16]  arrd  provide  simplified  expressions  as  a  fimction  of 
the  modified  Bessel  fmrctions  and  their  approximations.  Tire  applicatiorr  of  the 
corrhrruatiorr  method  to  the  somce  positiorr  estiuratiorr  problem  corrsists  of 
corrvolvirrg  the  origural  LS  objective  fmrctiorr  witlr  a  smootlrirrg  kenrel  before 
solvirrg  tire  locahzatiorr  problem.  Tire  idea  is  to  initially  apply  excessive  smootlrirrg 
to  the  exterrt  tlrat  the  objective  furrctiorr  becorrres  corrvex.  Murinrizmg  a  corrvex 
fmrctiorr  is  straightfor’wwd  to  inrplenrerrt  as  local  search  tecluriques  such  as 
gradierrt  descerrt  metlrods  become  global  search  methods.  By  progressively 
taper’urg  the  smootlriirg  effect  arrd  solvirrg  tire  problem  star  tirrg  from  the  solutiorr  of 
the  previous  step,  a  firral  more  robust  solutiorr  is  ofterr  obtairred.  If  suitable  care  is 
applied  to  the  taper’iirg  process  tlris  progressive  optmrizatiorr  method  usually 
results  irr  arr  accurate  approxinratiorr  of  the  R-LS  solutiorr. 

Tire  nraiir  issues  of  tlris  approach  are  the  choice  of  a  smootlrirrg  kenrel  arrd  tire 

der’ivatiorr  of  closed  form  expressiorrs  for  the  smootlred  objective  fmrctiorr  arrd  its 

2  -2 

gradierrt.  Usmg  the  Gaussiarr  kenrel,  g(ii,A)  =  e""  ,  Destirro  et  nl.  [16]  were  able  to 

der'ive  expressiorrs  for  the  smoothed  objective  furrctiorr  arrd  its  gradierrt  as 


arrd 


/=1 


/=1 


^3  r-  ^ 
2  A- 


N 


N 


v,/,(p) = Z'.-,v,/;(p) = 2;-v;(r„A) 


1=1 


/=1 

A 


P-P. 

I|P-P/II 


(26) 


(27) 


i2 

f 

■>  > 
3 

f .  •) 

5  rr 

e  ^ 

2,Fi 

-3,F, 

Tire  var’iable,  r.,  represerrt  the  i^h  serrsor-somce  Euclicharr  distarrce,  ||p-Pj||,  where 
p  is  the  mrkrrowrr  source  positiorr  hr  Tire  fmrctiorr  ^F^[a,b,z)  is  the  corrfluerrt 


hypergeometi’ic  furrctiorr.  Tire  paranreter,  A>0,  corrtrols  the  anromrt  of  smootlrirrg. 
Tire  larger  dris  paranreter  is,  the  rrrore  smootlrhrg  is  applied  to  the  objective 
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function.  As  X  approaches  zero,  tire  fimction  (26)  converges  to  the  original 
objective  function,  /(p). 


To  compute  the  smoothed  objective  hmchon  and  its  gradient,  we  depart  from  the 
original  derivations  given  in  [16]  arrd  take  advantage  of  two  useful  properties  of 
confluent  hr^rergeometiic  frmctions.  Tlrese  are 


3 

2 

(5 


-  \ 


2  - 


1 


2/2 


(2.-  +  3)/„(^)  +  (2.-  +  l)/.(^) 


(28) 


(29) 


2  2 

where  I^{.)  and  7j(.)  are  modified  Bessel  frmctions  of  the  first  kind  of  order  0  arrd 

1.  Usiirg  these  two  iderrfifies,  5'j(7',2)  reduces  to  5,(r,/l)  =  -e“"'(/„(.Y)  +  /,(.T))  where 

•) 
r' 


X  = 


2X- 


Now  to  compute  the  modified  Bessel  fmrcfiorrs,  it  is  advarrtageous  to  nrake  use  of 
pohmonrial  approxinratiorrs,  developed  by  E.  E.  Allerr  [22].  Tlrese  are 

4xe-^I„  (x)  =  0.39894228  +  0.0 1 32859 1 7 1  +  0.002253 1 87 1" 

-0.001575649t^  +0.0091628087'  -0.020577063f' 

+  0.026355372t®  -  0.016476329i"  +  0.003923767 =  Q^(x) 

^^xe-^I^{x)  =  0.39894228  -  0.039880242/  -  0.003620183r 

+  0.0016380147  -0.01031555/"  +0.022829673/^ 

-0.028953121/*  +0.017876535/"  -0.004200587/*  =^,(x) 

3  75 

where  /  = -  arrd  .y  is  fir  the  firter'r'^al  [3.75,  +  oo).  The  largest  approxinrafiorr  error 

.r 

fir  tlris  irrffirite  regiorr  is  =  2.0x1 0~".  Tlris  regiorr  is  of  firterest  wherr  A  approaches 
zero.  Tire  i“'  (mrweiglrted)  term  fir  the  smoothed  objective  furrctiorr  is  approxfiirated 
as 


/i(P)  =  4f‘  +iT  +X  '''((2.V,  +l)/o(Y)  +  2.r,/i(.v,)) 

Jlx, 


^  ^  ^  j  _ ^  j 

=  dr+ 1)- (1  +  — ) - (1  +  -—)Qo (.r, )  +  Qi (x, ) 


2.V, 


2y, 


(30a) 


where  .y,=— arrd  y,>3.75.  Its  coriesporrdfirg  gradierrt  becomes 

approximately 


2-V^^(^.(.y)  +  e>,(.y)) 


(P-P,) 


(30b) 
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For  smaller  values  of  x  (0<x<3.75),  which  correspond  to  large  values  of  X,  we 
make  use  of  the  following  approximations  [22] 

/„  (x)= 1  +  3.5 156229 +  3.0899424 + 1.2067492  i® 

+0.2659732i®+0.0360768i'°  +  0.0045813i'" 

x"'/,  (x)  =  0.5  +  0.87890594 +  0.5 1498869 +  0. 15084934  f 

+0.02658733i*  +  0.00301532i“’+0.0003241lP 

where  The  approximation  errors  are  less  than  » 3.0x10“**  for  /^Cx)  and 

» 1.0x10“**  for  /i(x).  In  this  region,  the  i*  (unweighted)  term  of  the  smoothed 
objective  function  becomes 

x;  (p)  =  Jf  +  r^{\  +  -^)  -  ((2x,  +  l)/„(x,)  +  2x,  A(x,))  (31a) 

2x,  ^2x, 

and  its  gradient  reduces  to 

Vp/;(p)=  2-^/^^7^e“"■(/„(x,.)  +  /l(x,.))  (p-p,.)  (31b) 

I  'i'  ) 

where  x,  =  and  0  <  x,  <  3.75  . 

'  22" 

Destino  et  al.  [16]  have  derived  an  expression  for  a  starting  initial  value  for  the 
smoothing  parameter,  2 ,  that  ensures  the  convexity  of  /^(p) .  This  value  is 

2  =^^max{d.}  (32) 

Starting  from  2^  and  an  arbitrary  initial  source  position  estimate,  p„,  a  practical 
gradient  descent  algorithm  would  converge  to  the  global  minimum  because  (p) 

is  convex  by  construction.  This  global  solution  becomes  the  starting  point  for  the 
next  gradient  descent  iteration  where  2  is  stepped  down  to  a  new  value.  By 
progressively  reducing  the  value  of  2  in  discrete  steps  and  solving  the  localization 
problem  starting  from  the  solution  of  the  previous  step,  a  succession  of  more 
accurate  solution  updates  are  usually  obtained.  The  main  aim  of  this  solution 
method  is  to  steer  the  final  solution  towards  the  global  minimum  of  the  original 
cost  function  /(p).  As  argued  in  [16]  this  is  almost  always  achieved  after  a  few 
fixed  iterations  of  2  .  In  our  computer  simulations  we  reduced  2  according  to  the 
geometric  progression  law,  2^^=  2^0.6*,  with  ^  =  0,1, 2,  •••,10  and  added  a  last 
iteration  where  2  =  0,  resulting  in  a  total  of  12  iterations. 
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6.  Geolocation  Using  the  Null  Space  Method 

Ill  tilts  section  we  focus  on  applying  Reid  and  Zlii's  method  [10]  (or  RZ  for  short)  to 
solve  systems  of  polynomial  equations  with  finitely  many  solutions  (i.e.  zero- 
dimensional  pohuomial  systems).  We  have  selected  to  use  RZ  because  of  its 
simplicity  and  efficiency  to  solve  tliis  problem.  Althougli  tlus  method  solves 
exactly  the  pohuomial  system  associated  witli  (1)  for  an  arbitraiy^  N ,  it  becomes 
computationally  prolnbitive  and  inaccmate  as  N  increases  beyond  3,  winch  coincides 
widi  die  minLmum  lumiber  of  measurements  rec|uh’ed  to  generate  an  unambiguous 
somce  position  solution  in  the  plane. 

Retelling  to  Figure  5  and  substituting  u.  =[cos^  ,  siii0.]^  in  (3)  and  expanding,  we 
obtain 

x  =  A  +  '^M)df  cosdf 

_  (33) 

v  =  B  +  y'  ii’d,  sill  6, 

•'  II  I 

N 

where  p  =  [.t  y]^ '  ^  =  X  ”  ^  “  X  • 


Figure  5:  Relationship  betioeai  p-p,  and  bearing  angle  6^ 

From  Figiue  5  it  also  follows  that  cos^  (y-y,)  =  sin^  (.r-.y),  winch  leads  to  N 

trigonometilc  equations  of  N  miknown  angles, 

a 

(5-  yy)cos0^  +(.Y^  -^)siii0^  siii(0, -^^.)  =  0,  (1  <  j  <N)  (34) 

1=1 

Wlien  the  number  of  sensors  is  restricted  to  N=3,  then  apphing  (34)  for  each 
beaiing  angle  leads  to  a  system  of  tluee  tilgonometric  equations  of  tlnee  variables. 
By  applying  the  nomiahzation  property,  c:  +  5/  =  1 ,  where  c,  and  s,  stand  for 
cos^  and  sin  6*  ,  we  obtain  die  equivalent  pohuomial  system 
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(S  -  )  Cj  +  (Xj  -  ^)  Sj  +  -  CjSj  )  +  w^d^  (SjCj  -  CjSj )  =  0,  cf  +  sf  - 1  =  0 

(S  -  )  Cj  +  (^2  -  ^)  Sj  +  Wj Jj  (SjCj  -  Cj^2 )  +  W3 J3  (S3C2  -  C3S2 )  =  0,  C2  +S2  -1  =  0  (35) 

(S  -  J3  )  C3  +  (X3  -  ^)  S3  +  Wj  Jj  (5jC3  -  CjS3  )  +  Wj  (■^2*^3  “  ^3^3  )  =  0,  C3  +  ^3^  -  1  =  0 

which  can  be  solved  using  a  number  of  polynomial  system  solving  techniques 
[8] [9] [10].  An  upper  bound  solution  count  for  this  system  is  2®  =64.  Using 
Groebner  based  analysis  the  solution  count  is  only  20. 

In  the  following  we  briefly  review  Reid  and  Zhi's  method  [10]  for  solving 
polynomial  systems  and  provide  an  algorithm  that  exactly  solves  (35),  subject  only 
to  numerical  round-off  errors.  A  polynomial  system,  S,  in  iR[Xj,X2,---,x^]  of  degree  q 

may  be  written  in  matrix  form  as  =0  where  =  [xf,xf■'x2,.••,^^,A’•••Amd]• 

ln  our  localization  problem  S  consists  of  the  six  polynomials  given  in  (35).  The 
number  of  variables  is  m  =  6  and  the  system's  degree  is  q  -  2.  The  main  idea 
behind  the  RZ  method  is  to  expand  the  polynomial  set,  S,  by  multiplying  each 
polynomial  of  S  by  monomials  up  to  some  total  degree,  which  is  determined  by 
the  involutive  criterion  for  zero-dimensional  polynomial  systems  [10].  The  process 
of  expanding  S  is  achieved  through  prolongation  [10].  A  single  prolongation  of  S 
increases  the  system's  largest  total  degree  by  one.  Prolongations  up  to  an  order  k 
result  in  the  polynomial  system,  which  can  be  represented  in  matrix  form  as, 
M,x,  =0,  where  the  non-zero  entries  of  M,  consist  of  shifted  polynomial 
coefficients  of  S  and  the  total  degree,  t,  is  related  to  prolongation  order,  k,  through 
t-q  +  k.  The  vector,  x, ,  stacks  all  monomials  up  to  total  degree,  t.  Clearly  all 
solutions  reside  in  the  null  space  of  M,  for  all  t>q.  The  null  space  is  represented 
in  matrix  form  as  B,  where  the  columns  constitute  a  null  space  basis  of  M, . 

In  addition  to  the  prolongation  operator  Reid  and  Zhi  defined  the  projection 
operator  which  is  applied  to  null  space  bases.  A  single  projection  involves 
removing  rows  of  B,  associated  with  the  total  degree,  t.  This  projection  step, 

denoted  as  ;^(S**^^)  has  the  effect  of  reducing  the  total  degree  of  x,  by  one,  hence 
mapping  x,  into  x,_j.  Higher  order  projections  are  represented  as  ;r^(S^*^)  where  £ 
is  the  projection  order  (0  <  ^  <  A: ).  The  main  result  of  the  involutive  criterion  is  that 
one  seeks  the  smallest  prolongation  order,  k,  such  that  there  exists  a  projection 
order,  £ ,  satisfying 


dim;r^(S**^^)  =  dim;r^'''(S^*^^)  =  dim  (36) 

where  dim  denotes  space  dimension.  If  for  a  given  order,  k,  more  than  one  £ 
satisfies  this  relationship,  then  we  select  the  largest  projection  order  £  [10].  The 
common  dimension  corresponds  to  the  number  of  solutions  of  S . 
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Table  I:  Reid  and  Zhi's  Algorithm  (Null  Space  Method)  [10] 

1.  Construct  the  expanded  sparse  matrix,  M,,  where  t  corresponds  to  the  involution 
criterion  degree. 

2.  Compute  the  null  space  of  M,  and  stack  its  basis  along  the  columns  of  B, . 

3.  Perform  I  and  f +  1  projections  on  B,.  In  other  words,  keep  rows  associated  with 
monomials  up  to  total  degree  t-£  and  t-£-l  and  call  the  resulting  matrices  B  and  Bj, 
respectively. 

4.  Compute  the  SVD  of  Bi  =  [Ui  U2]-[i:i  0]^-V"' 

5.  Form  the  multiplication  matrix,  ,  with  respect  to  a  selected  variable,  x,  ,  using 

=  U[  -B^.  -V-i:-'  where  B^  consists  of  rows  of  B  corresponding  to  the  monomials 
■ 

6.  Compute  the  eigenvectors,  v  . ,  of  and  recover  all  solutions  of  S  from  the  entries  of 
UjV^.  in  (after  proper  scaling) 

7.  If  the  polynomial  system,  S,  is  the  result  of  a  real  polynomial  optimization,  then  select 
the  globally  optimal  real  solutions 


Table  II:  Dimension  of  k' 


N=3 

k=2  (f=4) 

k=3  (t  =5) 

fc=4  (f  =6) 

k=5  (t  =7) 

e=o 

57 

63 

76 

92 

e=i 

26 

21 

20 

20 

e=2 

16 

21 

20 

20 

e=3 

6 

16 

20 

20 

e=4 

1 

6 

16 

20 

e=5 

0 

1 

6 

16 

e=6 

0 

0 

1 

6 

1=7 

0 

0 

0 

1 

Once  the  prolongation  and  projection  orders  associated  with  involution  are  known, 
Reid  and  Zhi  have  shown  how  to  construct  a  multiplication  matrix,  ,  with 

respect  to  a  given  variable,  x^,  and  apply  eigen-decomposition  to  determine  the 
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solutions  of  S .  It  turns  out  that  the  multiplication  matrix  can  be  generated  from 
knowledge  of  .  The  algorithm  to  compute  this  matrix  is  detailed  in  Table  I. 

Table  II  shows  the  dimension  of  for  our  polynomial  system  (35),  for 

different  prolongation  and  projection  orders.  The  first  row  gives  the  dimension  of 
the  null  space  of  the  expanded  matrix,  M, .  The  involutive  condition  (36)  is 
satisfied  at  t  =  6  {k  =  A)  and  f  =  2  as  indicated  by  the  boldface  entries  in  the  table. 
The  number  of  solutions  is  20.  The  size  of  the  expanded  matrix,  M,^g,  is  1260x924, 
which  is  rather  large  to  efficiently  compute  the  solutions  of  (35).  If  we  instead 
consider  column  t  =  5  in  the  table,  we  observe  that  dim ;r(S*^^)  =  dim ;r^(S^^^)  =  21 . 
The  associated  expanded  matrix  is  of  much  reduced  size  (504x462),  allowing 
a  much  faster  application  of  the  algorithm. 

Running  the  algorithm  on  a  PC  equipped  with  Intel(R)  Core(TM)i7-2600  CPU  @ 
3.40  GHz,  takes  about  0.03s  to  find  the  global  minimum  of  (1).  For  a  larger  number 
of  range  measurements  (TV >  4),  solving  (34)  is  in  principle  possible,  but  proves  to 
be  impractical  using  the  RZ  method  because  it  scales  poorly  with  the  number  of 
measurements.  Other  polynomial  solution  methods  also  suffer  for  larger  N. 
However,  the  computational  complexity  of  some  grows  much  faster  than  others. 
As  an  example  we  have  applied  the  polyhedral  homotopy  method  [8]  to  the  N=3 
and  N=4  localization  problems,  and  it  took  0.09s  and  0.5s  to  respectively  solve 
them  (20  solutions  for  N=3  and  64  solutions  for  N=A).  In  comparison,  the  RZ 
method  resulted,  on  average,  in  0.03s  and  13s,  respectively. 


7.  Computer  Simulations 

This  section  provides  a  comparative  performance  analysis  in  the  plane  in  terms  of 
source  localization  accuracy  and  processing  speed  for  the  position  estimators 
presented  in  this  paper.  The  R-LS  estimator  is  implemented  using  exhaustive  grid 
search.  The  purpose  of  this  analysis  is  to  distinguish  those  estimators  having 
desirable  properties  such  as  low  computational  effort  and  high  localization 
accuracy. 

All  computer  simulations  are  conducted  in  Matlab  using  a  total  of  5000  realizations 
for  each  experiment.  Each  experiment  consists  of  fixing  the  range  measurement 
standard  deviation,  a,  and  randomly  selecting  sensor  and  source  positions  within  a 
square  of  20  by  20  length  units.  This  allows  us  to  have  a  fair  comparison  between 
the  algorithms  regardless  of  geometry.  The  range  measurements  are  computed  as 
true  Euclidean  distances  plus  zero-mean  Gaussian  errors  with  standard  deviation. 
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a.  This  setup  enables  us  to  compute  an  average  mean  square  error  (MSB)  for  each 
estimator  across  the  same  5000  sensor/ source  positions  and  noise  realizations. 

Table  1  provides  MSB  results  for  the  simulated  source  localization  estimators  as  a 
function  of  the  measurement  standard  deviation  error,  a.  Bour  values  for  a  were 
chosen  as  shown  in  the  table.  The  number  of  sensors  is  fixed  to  5.  As  is  clear  from 
the  table,  the  MSB  of  SR-LS  is  consistently  larger  than  the  MSB  of  the  remaining 
estimators.  The  R-LS  estimator  performs  best  for  small  measurement  errors.  Bor 
larger  measurement  errors,  the  CPnBSP-SDR  accuracy  is  clearly  superior  to  those 
of  all  estimators,  including  R-LS.  The  range  weighted  SR-LS  is  of  particular  interest 
as  it  displays  acceptable  localization  accuracy  yet  low  computational  effort. 


Table  1:  Mean  square  error  as  a  function  of  range  measurement  error,  a,  for  N=5.  Also  included  are 
indicative  processing  times  in  milliseconds  using  a  PC  with  Intel(R)  Core(TM)2  Duo  @ 
2.53GHz  and  4GB  of  RAM 


N=5 

SR-LS 

(2ms) 

RW- 

SRLS 

(2ms) 

SDR 

(90ms) 

SLCP 

(190ms) 

CPnBSP 

SDR 

(70ms) 

Contin. 

(135ms) 

R-LS 

a=0.01 

2.0e-4 

1.4e-4 

2.6e-4 

1.6e-4 

1.8e-4 

1.4e-4 

1.4e-4 

a=0.1 

0.0199 

0.0137 

0.0308 

0.0170 

0.0177 

0.0501 

0.0137 

a=1.0 

2.7874 

2.0547 

2.4247 

2.1213 

2.0158 

2.1574 

2.1117 

a=2.0 

11.577 

10.198 

9.3889 

9.9240 

9.1428 

10.126 

9.7469 

To  have  an  appreciation  of  how  nonlinear  localization  approaches  compare  with 
linear  ones,  we  have  also  implemented  the  BLUB  estimator  (Best  linear  unbiased 
estimator  [4] [23])  whose  average  MSB  results  are:  5.6e-4,  0.051,  5.32  and  18.98  for 
a=0.01,  0.1,  1.0  and  2.0,  respectively.  Despite  its  simplicity  and  very  low  average 
processing  time  ('®0.23ms),  it  is  clearly  a  poor  localization  estimator.  Indicative 
processing  time  averages  are  also  included  in  Table  1  to  highlight  computational 
effort.  Notice  that  the  CPnBSP  SDR  exhibits  the  lowest  computational  effort  among 
all  relaxation  algorithms. 


Table  2:  Mean  square  error  as  a  function  of  range  measurement  error,  a,  for  N=10 


N=10 

SR-LS 

RW- 

SRLS 

SDR 

SLCP 

CPnBSP 

SDR 

Contin. 

R-LS 

a=0.01 

8.2e-5 

5.1e-5 

6.3e-5 

5.2e-5 

5.4e-5 

5.1e-5 

5.1e-5 

Q 

II 

o 

U 

0.0081 

0.0052 

0.0063 

0.0053 

0.0054 

0.0052 

0.0052 

a=1.0 

0.8108 

0.5916 

0.6361 

0.5521 

0.5551 

0.5413 

0.5413 

a=2.0 

3.3730 

2.8998 

2.6818 

2.4003 

2.3759 

2.3404 

2.3446 
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Table  2  provides  MSB  results  for  N=10.  The  SR-LS  provides  the  worst  localization 
accuracy  performance.  The  R-LS  and  the  continuation  method  appear  to  be  the 
most  accurate  estimators.  The  SLCP  and  CPnBSP  SDR  accuracy  performances 
come  next  and  are  close  approximations  to  R-LS.  The  CPnBSP  SDR's 
computational  efficiency  is  moderate  ('®95ms).  The  weighted  SR-LS  offers  a  slightly 
degraded  accuracy  but  a  much  lower  computation  time  ('»2.1ms),  making  it  an 
attractive  estimator  for  practical  use. 

The  next  set  of  computer  simulations  aim  at  comparing  source  localization 
performances  for  fixed  sensor  and  emitter  positions.  Table  3  displays  the  MSB  of 
six  estimators  using  a  fixed  scenario  of  5  sensors  positioned  at  (6,  4),  (0,  -8),  (-7,  3), 
(8,  -5),  (3,  -3).  The  true  source  location  is  maintained  at  (-2,  3).  The  MSB  associated 
with  the  Cramer-Rao  Lower  Bound  (CRLB)  is  depicted  in  the  last  column. 

Table  3:  Mean  square  error  as  a  function  of  range  measurement  error,  a,  for  a  fixed  scenario  of  5 
sensors.  The  last  column  shows  the  lowest  achievable  MSB 


N=5 

SR-LS 

RW- 

SRLS 

SLCP 

CPnBSP 

SDR 

Contin. 

R-LS 

CRLB 

a=0.01 

1.2e-4 

l.Oe-4 

l.Oe-4 

l.Oe-4 

l.Oe-4 

l.Oe-4 

l.Oe-4 

a=0.1 

0.0125 

0.0102 

0.0104 

0.0102 

0.0102 

0.0102 

0.0102 

a=1.0 

1.2697 

1.0769 

1.0654 

1.0563 

1.0512 

1.0512 

1.0212 

a=2.0 

5.1984 

4.7314 

4.4505 

4.4720 

4.6983 

4.6983 

4.0847 

As  is  clear  from  Table  3,  the  R-LS  performance  approaches  the  CRLB  for  relatively 
small  distance  measurement  errors.  It  clearly  departs  from  CRLB  at  a=2.0. 


Table  4:  Mean  square  error  as  a  function  of  range  measurement  error,  ct,  for  N=3 

MSB  Relative  MSB 


Range 

error 

R-LS 

CPnBSP 

SDR 

SDR  [7] 

CPnBSP 

SDR 

SDR[7] 

a=0.01 

0.176 

0.192 

0.21 

0.027 

0.05 

Q 

II 

o 

2.152 

2.103 

2.14 

0.27 

0.515 

a=1.0 

23.84 

21.78 

21.03 

1.373 

3.019 

a=2.0 

45.83 

43.11 

41.16 

2.016 

4.343 

Table  4  provides  MSB  results  when  N  is  restricted  to  3.  Bour  values  for  a  were 
chosen  as  shown  in  the  table.  The  R-LS  position  estimator  is  computed  using  RZ 
method  as  detailed  in  section  6.  Bor  very  small  range  measurement  errors,  the  R- 
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LS  estimator  performs  best.  For  large  measurement  errors,  the  accuracy  of  the  SDR 
algorithm  [7]  is  marginally  superior  to  that  of  CPnBSP  SDR.  The  last  two  columns 
in  the  table  show  relative  MSB  with  respect  to  R-LS  estimates  (i.e., 

}  and  -P^j-r^ir) )  and  clearly  highlight  the  superior 

approximation  accuracy  of  the  CPnBSP  SDR  to  R-LS. 

To  summarize,  there  is  not  much  to  gain  from  the  SLCP  estimator  as  its  MSB 
performance  is  comparable  to  that  of  CPnBSP  SDR  but  its  computational  effort  is 
larger.  The  SDR  estimator  is  generally  less  accurate  than  CPnBSP  SDR.  The  SR-LS 
is  the  least  accurate  estimator.  The  continuation  method  can  be  useful  but  its 
performance  is  less  reliable.  Its  major  drawback  is  that  it  does  not  provide  a 
convergence  certificate  to  the  global  minimum.  Grid  search  R-LS  is  not  practical 
and  does  not  guarantee  convergence  to  the  global  minimum.  The  RZ  method 
provides  an  exact  R-LS  solution  and  is  reasonably  fast  for  small  localization 
problems  where  N  =  3.  To  conclude  it  appears  that  both  the  CPnBSP  SDR  and  the 
weighted  range  SR-LS  are  the  preferred  estimators  for  arbitrary,  N.  Both  may  also 
be  used  for  source  localization  in 9?^.  Bor  time  critical  applications  the  range 
weighted  SR-LS  is  obviously  the  position  estimator  of  choice  as  it  offers  the  best 
compromise  between  estimation  accuracy  and  computational  effort. 


8.  Conclusions 


This  report  presents  a  review  of  range-only  localization  algorithms  that  aim  to 
approximate  the  asymptotically  efficient  but  computationally  demanding 
maximum  likelihood  estimator,  which  reduces  to  R-LS  if  the  range  measurement 
errors  are  normally  distributed. 

The  optimization  criteria  employed  by  the  range-only  localization  algorithms  cover 
semidefinite  relaxation,  global  continuation,  computer  algebra  and  the  generalized 
trust  region  subproblem  method.  Bxtensive  computer  simulations  were 
implemented  to  compare  not  only  the  estimation  performance  of  the  different 
algorithms,  but  also  their  computational  effort.  The  new  proposed  algorithms 
presented  in  this  paper  were  observed  to  perform  well.  The  RW-SRLS  estimator  is 
suitable  for  real-time  applications.  Its  accuracy  is  only  slightly  worse  than  that  of 
R-LS.  The  CPnBSP  SDR  algorithm,  comes  second  in  terms  of  computational  effort 
but  achieves  an  improved  MSB  performance,  which  usually  matches  that  of  R-LS. 
The  CPnBSP  SDR  seems  to  slightly  outperform  R-LS  in  scenarios  where  only  a 
small  number  of  measurements  are  available  or  the  measurement  errors  are  large. 
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